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Abstract - We calculate the linear evolution for a class of mixed dark matter models, where the hot 
component derives from the decay of a heavier particle (CntHDM models). These models differ from 
ordinary mixed models based on massive neutrinos (CHDM models), for which the hot component has a 
phase space distribution which derives from a thermal one. In CntHDM models the density of the hot 
component and the derelativisation redshift of its quanta are indipendent parameters. In this work we 
provide the spectra for a number of CntHDM models, and compare them with CDM and CHDM spectra. 
If PQ symmetry and SUSY simultaneously hold, the lightest standard neutralino can be expected to decay 
into axino and photon. We briefly summarise the features of this particle model which gives rise to the 
cosmological framework discussed here, as a fairly generic consequence. Other frameworks which lead to 
similar models are also briefly discussed. 
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1. Introduction. 

Mixed dark matter (MDM) models (Valdarnini & Bonometto 1985, Bonometto & Valdarnini 1985, Achilli 
et al. 1985, Holtzman 1989, see also Fang, Li, & Xiang 1984) allow a better fit of large scale data 
than most other models (Schaefer, Shah, & Stecker 1989; van Dalen & Schaefer 1992; Schaefer & Shafi 
1992; Taylor & Rowan-Robinson 1992; Davis, Summers, & Schlegel 1992, Holtzman & Primack 1993; 
Pogosyan & Starobinsky 1993; Liddle & Lyth 1993, Klypin et al. 1993, Klypin, Nolthcnius & Primack 
1993, Bonometto et al. 1994). Recent treatments of MDM were mostly based on a mixture of CDM, 
HDM and baryons with Sl c /r2/j,/Slfe = 0.6/0.3/0.1 (O: density parameters), although different ratios were 
also considered. If HDM is made of fermions of mass rrih, with gf spin states, which decoupled when 
T » rrihc 2 , it is Qh — 0.3(ra/ig//14cV). The standard candidate as HDM component is the r-neutrino, 
assumed to have a mass rrih = 7 eV. 

Hereafter this version of MDM will be called CHDM. Within its context tth and the derelativization 
rcdshift for the hot component Zd e — m^c 2 /3T . U are connected through rrih (T , y is the present neutrino 
temperature). In turn, Zd e fixes the mass scale M^^ above which fluctuations in the hot component are 
not erased by free streaming after their entry in the horizon. 

Instead of being composed of massive neutrinos, the hot component could have arisen from the decay 
of heavier particles. In this case, the link among rrih, gf, &h, za e , and Mp^ is not so cogent. There are 
still precise constraints due, e.g., to primeval nucleosynthesis, but the parameter space is spanned by two 
independent continuous parameters (e.g., and Zde), instead of a continuous and a discrete parameter 
(e.g., and gf). The distribution of hot particles in momentum space is then quite different from the one 
originating from a thermal distribution. In this note we shall discuss such alternative to CHDM (hereafter 
denominated CntHDM). 

The decay of mother particles takes place when fluctuations over the relevant scales are still well 
outside the horizon. Daugther particles are then relativistic and an analytical treatment of the evolution 
of fluctuations during these stages can be given. The critical events to set up the shape of the final 
spectrum are the entry of the fluctuations in the horizon and the derelativization of the particles forming 
its hot component. These two events can take place in different order, according to the scale considered 
and to the model parameters. At a redshift ~ 10 3 , primeval plasma (re)combines. In purely baryonic 
models recombination is a critical step, as it sets up a minimal mass scale for the survival of baryonic 
fluctuations and provides information on very-small-scale CBR fluctuations. In models dominated by DM, 
these stages are not so essential, as baryons will later fall in potential wells created by CDM and/or HDM, 
while the relevant scales for CBR fluctuation observations however enter the horizon after recombination. 
The presence of radiation and baryons is critical for fluctuations entering the horizon before equivalence. 
The growing of the CDM component is then almost stopped by Meszaros effect, i.e. by the fact that 
baryon and radiation can only oscillate as sound waves while, until equivalence, they are the dominant 
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substance in the Universe. After equivalence, as is known, DM fluctuation growing can restart. 

Furthermore, the baryonic component has a non neglegible influence on the detailed final spectrum. 
Our models take also into account a massless neutrino background and its fluctuations (which are rapidly 
dissipated at the entry in the horizon). Three, two or one massless neutrinos were considered, according 
to the model parameters, as is discussed below. 

The main output of this paper are a number of spectra for CntHDM models, computed in detail 
through all their linear evolution, down to z = 0. The models considered are given in Table 1 of sec. 2. As 
is also outlined in sec. 4, the study of each model still requires a fairly large amount of CPU, and further 
numerical analysis is still in progress. 

A particle model, recently discussed by Bonometto, Gabbiani & Masiero (1994, BGM), originates 
CntHDM models as a fairly generic consequence. We shall briefly outline some features of this particle 
model in sec. 2. The class of cosmological models described here is however relevant, also independently 
of that specific particle physics. In sec. 3 we will study the analytical part of the evolution of density 
fluctuations. In sec. 4 the numerical treatment will be described. Sec. 5 is devoted to a discussion of the 
results. 

2. A particle model. 

In this section we debate a specific particle model which gives rise to a cosmological framework of the 
kind discussed in the next sections. Although this model is complitcly consistent, the cosmological results 
have a validity which extends well beyond this peculiar particle physics. In recent years cosmologists have 
become familiar with supersymmetry (SUSY) and Peccei-Quinn (PQ) symmetry. Both of them provide 
CDM candidates. The lightest mass eigenstate of neutral fermion SUSY partners, called neutralino (\), 
has a mass m x ^ 30 GcV and must therefore decouple when T is significantly below m x (Lee & Weinberg 
1977, Ellis ct al., 1984). The goldstone boson associated to PQ symmetry, arising at T < Vpq (scale of 
PQ symmetry breakdown), is called axion (a). If both SUSY and PQ symmetry are implemented, also 
SUSY partners of a's exist; the fermion partner of a is called axino (a). In BGM it is shown that we can 
expect two a background to exist. Besides thermal a's, decoupling slightly below Vpq, probably diluted 
by entropy inflow(s) at later phase transition(s), non-thermal a's originate from \ decay. The former 
component is cold, the latter one is a non-thermal hot component. It is possible that all DM consists of 
a's, which can account for both components of MDM. 

The decay \ ~~ > a + 7 takes place at z x ^ y ~ 10 9 and a's have initially momentum P = m x c/2. They 
derelativize at Zd e ~ z x ^ y 2ma/m x . The outcoming density parameter fla — flh is linked to the axino 
mass (rria) and \ number density at the eve of their decay (n x ). These quantities are linked by a number 
of particle relations, discussed in BGM and summarized in fig. 1. In such figure m s f is the scale of the 
soft SUSY breaking (sfermion scale. 

Different curves refer to possible choices of m x and m t (top quark mass), as is detailed in its caption. 
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Table 1: Parameters for CntHDM models 



model 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 




0.30 


0.30 


0.30 


0.15 


0.20 


0.15 


0.20 


0.20 


0.15 


0.20 


Wio 4 


1.5 


6 


3 


1.5 


0.7 


0.7 


0.9 


1 


0.4 


0.4 


massless z/s 


2 


3 


3 


3 


2 


2 


2 


1 


2 


1 



For each value of m 8 f, fig. 1 provides the values allowed for the axino mass ma (decreasing curves) and 
for the quantity ZdeVpQ/Qa, where the particle parameter Vpq, expected to range between 5xl0 9 GeV 
and 10 12 GeV, is combined with Zd e (a derelativization redshift) and fij (hot DM density parameter). 

Fig. 1 shows that the quantity z de VpQ/fla can be expected to range between 10 13 GeV and 10 16 GeV, 
for reasonable values of m s f. In turn, it is possible to find various combinations of particle parameters 
which give Oa in the range 0.1-1, and Zd e in the range 10 3 -10 5 . 

Further constraints are set by the number of massless neutrinos still present during primeval nucleosyn- 
thesis. Let g v be the number of massless v spin states present during the nucleosynthesis (z/s with mass 
much below nucleosynthesis temperatures are to be considered massless). Standard stellar abundance 
data for light elements allow up to g v ~ 7. Accordingly, if g v = 6 (3 standard neutrinos) we are allowed 
g + = 1 extra spin states. It can be g + = 3 or g + = 5, if one or two neutrino(s) are supposed to be quite 
heavy and to have already decayed before nucleosynthesis. Neutrinos with a mass m v =7 cV (= 8 x 10 4 K) 
derelativize at a redshift z ~ 1.5 x 10 4 and yield fi„ = 0.3. A particle allowed half of their energy density 
at nucleosynthesis and derelativizing at Zde will then yield fi/, ~ Zde/10 5 . If the number of extra neutrino 
spin states allowed during nucleosynthesis is g + , the limit for flh reads: 

fl h < g+ z de /10 5 . (2.1) 

The number g v of effective massless neutrinos is used in the equations ruling the expansion of the Universe 
through the quantity w v = (7/4)(4/ll)*g v , that will be used below. 

For all our models we took Sltotal = 1, fi& = 0.05, and H = 50 kms _1 Mpc _1 . Values of fi/, and 
Zde ranging in the intervals 0.15-0.30 and 4 x 10 3 -6 x 10 4 were considered. Models with low z de show 
interesting features and the parameter space was explored in that region in more detail. In Table 1 we 
report the values of fi/,, Zd e , and the number of massless neutrinos, for those models whose spectra are 
worked out here. 

3. Density fluctuation evolution: analytical aspects. 

Let f(x a ,Pa,t) be the ntH (non-thermal hot) a distribution in phase space (x a , P a = Pn a : particle 
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coordinates and momentum). We must also consider the source of a's due to x _ decay: S(x a , P a , t). Up 
to first order, 

f = f (P,t)[l + e(x a ,P,t)], S = S (P,t)[l + a s (x a ,P,t)] (3.1) 
(Notice that I and o s do not depend on n a ). The analysis is based on the equation 

s = df/dt + (df/dx a )x a + (df/dP)P (3.2) 

taking into account that x a = cn a {P/ P )(T /T), P = P[(T/T) + y/A\. Here P Q c is energy, T is radiation 
temperature, y yields perturbation self-gravity and will be given in eq.(3.7) below. At the 0-th order, 
cq. (3.2) yields 

S = df /dt + P(f/T)(df /dP) , (3.3) 
while the source term, due to x _ decay, is 

S = (h 3 /4Tr)n x , dg {T/T x . Dy ) 3 exp(-t/t XiDy ) {t x , Dy P 2 )- l 8{m x , Dy c/2 - P) (3.4) 

{tx,Dy'- X mean life, n Xtdg :xs number density at the decoupling ). With such source terms the solution of 
cq. (3.3) reads 

f a - (Y?j^)n xM exp(-Q) (2Q/P 3 ) 6{m XiDy c/2 - P) (3.5) 

with 

Q = (2P/m x cf(t/t x . Dy ) (3.6) 

(notice that Q S(m x c/2 — P) = t/t x ,D y )', X decays occur in radiation dominated era, when t/t x ,jj y — 
{T x . Dy /T) 2 . 

We shall treat the first order equation for the case of a plane wave in x% direction. In such direction 
the component of P is /i P. Let also h a p — ri a p — g a p express the deviation of the metric tensor g a p from 
purely spatially flat r\ a p components (a, /3 = 1..3; t being the universal time, only g ao does not vanish, 
besides of g a (})- Then 

y = (1 - fi 2 )3h t - (1 - 3fi 2 )h 33 = 2P (fi)h t + 2P 2 (fi)h 3 , (3.7) 

where 

3 

3h t = ^2 h aa , h 3 = h t - h 33 . (3.8) 

a=l 

Let us then take e = e(k, P, fi, t) exp(ikx 3 T/T a ), in order to treat each length scale separately (p c : present 
critical density; Pi: Legcndre Polynomials). 

The first order equation splits in two parts. Let be K — kcT/T a . For P < m x c/2, 

e + iefiKP/P = (1 + 2Q)y/4 , (3.9) 
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while, for P = m x c/2, it is e = a s + yt/2. The latter equation holds at x decay and sets initial conditions 
to eq. (3.9). Such equation is to be solved for any \x and P; however P is related to Q (see 3.6) and the 
distribution f Q decreases exponentially with Q (see 3.5). This will allow to select a finite set of Q values. 
As far as p, is concerned, we perform the expansion 

oo 

e(k,P, f i,t)=J2(-i) l °i(k,P,t)Pi(ri . (3.10) 

Then, taking eq. (3.7) into account, eq. (3.9) yields the system 

o = -(l/3)tf(P/P>i + (1/2 + Q)h t , <7! = K(P/P )(a - 2a 2 /5) 

6-2 = K(P/P )(2a 1 /3 - 3a 3 /7) + (1/2 + Q)h 3 , a t = tf(P/P )(i_<7,_i - l+a l+1 ) (I > 2) (3.11) 

[here I- = 1/(21 - 1), 1+ = (I + 1)/(2Z + 3)]. 

Radiation is to be treated very much alike. The critical difference is that photons are masslcss and 
therefore all P behave in the same way (furthermore, of course, the P distribution is different). Integrating 
the photon distribution over P and expanding the outcoming S(k, fi, t) — X)/^o( — t)Pl{^)^ we obtain 
a closed system analogous to (3.11), which reads 

S = -{1/3)KS 1 + 2h t , S x = K(S - 2S 2 /5) - (ft - Aw)v 

S 2 = K{2S 1 /3 -3S 3 /7) + 2h 3 -S 2 9i//10 , Si = K{l_Si-t - l+S l+1 ) - S t v {I > 2). (3.12) 

Apart of the replacement 1/2 + Q — > 2, the main difference are the terms accounting for matter-radiation 
collisions: v = n e (TTC is the inverse 7 collision time (n e : electron number density, ar- Thomson cross- 
section) and w is the velocity field in baryons. Eqts. (3.11) and (3.12) are to be implemented by the 
equations for baryon and CDM density fluctuations (5b, S c ), metric perturbations (h t ) and t dependence 
of T. Such equations read: 

(f/Tf = T 3 [n b + (1 + w n u)T/T e + n h e h + fl c ] , 

S b = -Kw + 3h t /2 , w = wf/T + (S 1 - Aw)(p r / ' p b )v /3 , 

h = 2h t f/T + T 3 [S b n b + 2(S + S Q}V )T/T e + e Q h + 8 C Q C } , 

h 3 /2 - -ht + (T 2 /kc)[wQ b + (Si + S ljV )T/3T e + ei Q h } . (3.13) 

Here p r /p b is the radiation/baryon density ratio, Sl c T e — 10 4 Q c T o is the temperature when CDM and 
radiation have equal densities, T n = (8ir/3)Gp c (T/T ) n , and 




(3.14) 
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Here % = (T Xi d v /T) 2 and, in most cases, it cannot be numerically distinguished from oo. X^e is the typical 
derelativization temperature of axinos. 

We took into account also massless neutrinos. Their fluctuations are rapidly damped as soon as they 
enter the horizon. In principle, to follow such damping in detail, the same number of equations as for 
radiation are to be used. Some approximations, introduced when residual neutrino fluctuations are below 
the precision limit of the whole algorithm, have been set up in order to limit the number of harmonic 
components. They will not be discussed here. In eq.(3.13), 8 , v and Si. u indicate the first two harmonics 
for neutrinos, multiplied by a factor which depends on the number of massless neutrinos in the model. 

4. Density fluctuation evolution: numerical aspects. 

We performed a set of numerical integrations of the integro-differential system of equation (3.11), (3.12), 
(3.13), plus massless neutrinos. The integral part of the system was treated by taking 10 values of 
P associated with those values of Q needed to perform the integrations in eq. (3.14) with the Gauss- 
Laguerre projection procedure. The number of p in eqts. (3.11)-(3.12) was self regulated to obtain a 
precision of 1:10 5 . The algorithm allowed up to a maximum of 499 harmonics, which were never reached. 
The integration routine is a distant descendant of the one used by Bonometto, Caldara, & Lucchin (1983), 
although the treatment of collisionlcss particles with ntH spectrum required substantial changes. 

Altogether the differential equations to integrate numerically are up to 5505. For each model 9 lenth- 
scales were considered (1 Mpc, 5 Mpc, 10 Mpc, 20 Mpc, 50 Mpc, 100 Mpc, 200 Mpc, 500 Mpc, 1000 
Mpc), approximately spanning the mass-scale range from 10 12 Mq to 10 21 Mq. Each model requires 
approximately 37 hours of CPU at a 4000/90 vaxstation. Small length-scale cases are obviously taking 
much more time than great length-scale ones. 

For the sake of example, in fig. 2 we report the detailed behaviour of density fluctuations for the 
various components, for the model 5 and the 100 Mpc case. Some irregularities are present in the plot of 
matter and radiation sound-waves, due to a number of prints not fully adequate to follow the oscillations 
in detail. 

In fig. 3 and 4 we plot the outcoming spectra at z = 0, arising from primordial Zel'dovich spectra, for 
the models 2 and 6. In the same plots the spectra for a pure CDM and a standard CHDM model are 
also shown. The CDM spectrum has also been worked out with our algorithm by setting = 10~ 3 . The 
discrepancy between such spectrum and the one given by Holtzman (1989), for the same cosmological 
parameters, is typically ^ 0.5%, and nowhere exceeds 2%. 

In order to argue on the physical interpretation of the outcoming spectra, a direct comparison with 
CDM and CHDM is more suitable than the detailed spectral shape. Henceforth we preferred to plot the 
differences Alog[P(fc)] = log[P(fc)] — \og[P(k)]cDM,CHDM , for the 10 cases we evaluated, rather P(k) 
themselves. In fig. 4, 5 and 6 such A log[P(fc)] are shown for the models 1-5, 5-7 and 8-10, respectively. 
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Table 2: Transfer function coefficients for CntHDM models 



model 


Cl 


C2 


C3 


c 4 


1 


-0.12991D+01 


0.11015D+02 


0.34827D+02 


0.40777D+03 


2 


-0.45687D+01 


0.60267D+02 


-0.13429D+03 


0.37425D+03 


Q 
o 




fl ^87^911-1-09 






4 


-0.17309D+01 


0.17613D+02 


0.69582D+02 


0.13268D+03 


5 


-0.11754D+01 


0.49944D+01 


0.13379D+03 


0.16949D+03 


6 


-0.20006D+01 


0.16444D+02 


0.10363D+03 


0.11294D+03 


7 


-0.13807D+01 


0.10134D+02 


0.10300D+03 


0.18045D+03 


8 


-0.12392D+01 


0.10929D+02 


0.86743D+02 


0.16254D+03 


9 


-0.22635D+01 


0.20171D+02 


0.10958D+03 


0.11857D+03 


10 


-0.19229D+01 


0.11527D+02 


0.14706D+03 


0.15760D+03 



5. Discussion. 

The main outputs of this work are reported in the figures and in Table 2, where we provide the coefficients 
for a numerical interpolation of the transfer functions 

T(k) = yfpWk = Ml [\ + J2 cj [fc/CMpc- 1 )]^ 2 j (5.1) 

through the coefficients Cj (the normalization coefficient TV will not bew given here). 

In fig. 5, fig. 6 and fig. 7 we compare the spectra of the models studied here and those of CDM and 
standard CHDM. Most features reported on these figures can be qualitatively understood taking into 
account the corresponding values of flh and Zd e - 

Model 1 (in fig. 5, solid curve). The behaviour of this model is not far from standard CHDM. Its 
transfer function nowhere exceeds CHDM by more than a factor 1.4. No attempt was made to fine tune 
f2/j and Zd e , to obtain a still nearer spectrum, but a better fit is possible. 

Models 2 and 3 (in fig. 5: dashed curve and dotted curve). They are characterized by a greater Zd e 
and the same f2/j, in respect to CHDM. Their transfer functions are then greater than CHDM over all 
those scales which enter the horizon when the hot component is already derelativized. At large k, the 
curves fall down towards CHDM, as is to be expected. 

Model 4 (in fig. 5: short-dashed dotted curve). This model differs from model 1 for a smaller flh (zde 
is the same). It is therefore intermediate between CDM and CHDM. The spectrum bends at the same 
scale as for CHDM and the transfer function ratio is then decreasing fairly slowly, in respect to CDM. 
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Models 5-6-7 (in fig. 6). Their behaviours show that a decrease of Zd e can widely compensate an 
increase of flh- The spectrum falls more rapidly than model 4 and begins to flatten at a smaller k. These 
models were considered to explore the parameter space with intermediate values for Zd e - 

Models 8-9-10 (in Sg.7). While model 8 still reproduces a situation similar to CHDM, just with 
smaller density for the hot component, the other two models (with low Zd e ) show a consistent decrease of 
the spectrum if compared to CDM for small ks and a stabilization of their behaviours at k/Mpc~ x ~ 1. 
This is likely to cause a significant improvement in respect to CHDM for the scales over which objects 
can form at high z. 

The linear evolution of primeval fluctuations can be used to evaluate a number of observational quan- 
tities that will not be fully enumerated here. Each one of them can be obtained by performing suitable 
integrations over particular regions of the spectrum. After COBE evaluation of 5T/T, it can be conve- 
nient to compare different models by requiring them to give the same result over the COBE angular scale. 
This is approximately obtained here by requiring that all spectra coincide at 500 Mpc. A more detailed 
integration of the numerical outputs for radiation can be performed to improve 5T/T results; the radiative 
component is evaluated by our algorithm down toz~ 500, and we shall use all relevant harmonics to give 
the expected 5T/T, over a wide range of scales, in a forthcoming work. 

Among the advantages of CHDM in respect to CDM is the behaviour of the spectrum over the scales 
yielding bulk velocities. However, the faster decrease of CHDM also leads to predicting a rather unadequate 
amount of high - redshift objects. This can be a problem for high redshift galaxies and QSO's (see, e.g., 
Pogosyan e Starobinsky 1994). More recently, it has been pointed out that the problem can be even more 
severe for objects associated to damped Ly Q clouds (Mo and Miralda-Escude', 1994, Subramanian and 
Padmanabhan, 1994, Kauffmann and Chariot, 1994, Ma and Bertschinger, 1994, see however also Klypin 
et al., 1994). 

In this connection we wish to draw the attention on the characteristics of the models 5-6 and 9-10. 
Such choice of parameters give a decrease of P(k) (in respect to CDM) down to intermediate scales, while 
the difference between CDM and such models tends then to stabilize, instead of further decreasing as in 
standard CHDM. More detailed computations are needed to verify up to which degree this improving of 
outcoming spectra can ease the fit with observational data. 

There can be scarse doubts that, thanks to the presence of an extra parameter, CntHDM models are 
more flexible than CHDM ones. However such extra parameter is connected with precise data of particle 
physics, which can be tested through experiments which are either already feasible or can be expected to 
be performed in a fairly nearby future. This does not only concern the axino model discussed in sec. 2, 
but also other models with decaying particles. Among them, models in which there can be decays of one 
neutrino flavour into lighter neutrinos can be of particular relevance. Although these decays have a more 
complex kinematics, the cosmological behaviour due to the outcoming spectra can be very similar to the 



9 



ones discussed here. 
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FIGURE CAPTIONS 



Fig. 1 - Particle and cosmological parameter for the a model. The four lines refer to different choices 
of m x and m top . The solid lines are for m x = 30GeV and m top — 120 GeV. The dotted lines are for 
m x = 30 GeV and mtop — 180 GeV. The short dashed lines are for m x — 60 GeV and rritop = 120 GeV. 
The long dashed lines are for m x = 60 GeV and m top = 180 GeV. The area between solid and long dashed 
lines is therefore however allowed for a combination of m x and rritop, hi their allowed ranges. 
Fig. 2 - Examples of time evolution of the different components (CDM, a's, massless z/s, baryons, 
radiation); t ri f = 1.44 x 10 n s. The ordinate scale is arbitrary. 

Fig. 3 - Spectrum at z=0 for the model 2 (see text). Dotted lines refer to a CDM and a CHDM model 

(the latter has Q. c /Q. h /fl b = 0.6/0.3/0.1). 

Fig. 4 - The same as fig. 3 for the model 6 (see text). 

Fig. 5 - The difference Alog[P(fc)] between the models 1 to 4 and CDM (lower curves) or CHDM (upper 
curves) is given in function of k. The solid curves refer to model 1. The dotted curves refer to model 2. 
The dashed curves refer to model 3. The short-dashed dotted curves refer to model 4. 
Fig. 6 - The same as fig. 5. Here the models 5 to 7 are considered. The solid curves refer to model 5. 
The dotted curves refer to model 6. The dashed curves refer to model 7. 

Fig. 7 - The same as fig. 5. Here the models 8 to 10 are considered. The solid curves refer to model 8. 
The dotted curves refer to model 9. The dashed curves refer to model 10. 
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